home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / ode_forward.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-17  |  6.9 KB  |  213 lines
  1. /*
  2. ### standard simple driver for integration ###
  3.  
  4. -----------------------------------------------------------------------------
  5. NOte:    Assume that choose_algorithm will handle nonsymp models
  6.     or models with odd dimensionality
  7.     -Integration is done in standard coods, that is, the coordinate system defining
  8.     the original vector field.
  9.     - Plotting is done in window coordinates, the coord system as defined in
  10.     window_size panel 
  11.     - Section is computed by section variables, the coordinates defined in section
  12.     panel
  13. -----------------------------------------------------------------------------
  14. - int_driver:    0 no sense of absolute time, integration from i=0 to i_max
  15.         1 knows absolute time, integration from int_start to int_end
  16.         in int_nsteps;
  17.         2 knows absolute time, integration from int_start to int_end
  18.         with variable steps but error tolerance less than int_eps 
  19. -----------------------------------------------------------------------------
  20. THIS DRIVER DEALS ONLY WITH CASE 0 & 1
  21. -----------------------------------------------------------------------------
  22. */
  23. #define TINY 1.e-10
  24.  
  25. int ode_forward()
  26. {
  27.     int i,j,iter,i2,dim2,choice,max_iter;
  28.     double fabs(),sin(),fmod(),new_time,new_time_step;
  29.     double *vx,*vx1,*vx2,*vx3,*var_sec,*var_sec_old,*vfull,*delx,*dvector();
  30.     void free_dvector();
  31.     extern int stop,nok,nbad,linear_interpolation;
  32.     extern int model,var_dim,full_dim,int_algorithm,nofrill,i_max,n_poincare,draw_all,polar_coord;
  33.     extern int polar_section,section_index,linear_interpolation, section_count;
  34.     extern int int_driver,int_nstep;
  35.     extern double int_start,int_end;
  36.     extern double cutoff,time,time_step,section_constant,average_return_time;
  37.     extern double *win_var_i,*win_var_f,*param;
  38.     extern int (*f_p)();
  39.     extern char string[];
  40.  
  41.  
  42.     stop = 0;
  43.     /* memory allocation */
  44.     dim2 = var_dim/2;
  45.     vx = dvector(0,var_dim-1);
  46.     vx1 = dvector(0,var_dim-1);
  47.     vx2 = dvector(0,var_dim-1);
  48.     vx3 = dvector(0,var_dim-1);
  49.     var_sec = dvector(0,full_dim-1);
  50.     var_sec_old = dvector(0,full_dim-1);
  51.     vfull = dvector(0,full_dim-1);
  52.     delx = dvector(0,var_dim/2);
  53.     
  54.     /* from window variables to primary variables */
  55.     from_window_variables(vx,win_var_i,polar_coord);
  56.     /* plot and initialize full window variables */
  57.     to_full_variables(vfull,vx,polar_coord);
  58.     if(draw_all){
  59.         (void) draw_record_orbit(vfull,0,0);
  60.     }
  61.     
  62.     /* initialize section variables and full variables including function variables */
  63.     if(!draw_all){
  64.         to_full_variables(var_sec,vx,polar_section);
  65.         for(i=0;i<full_dim;i++) var_sec_old[i] = var_sec[i];
  66.     }
  67.     
  68.     /* choose an appropriate algorithm */
  69.     choice = (int) choose_algorithm();
  70.     
  71.     /* get maximum number of iteration depending on the driver */
  72.     if(int_driver == 0)
  73.         max_iter = i_max;
  74.     else if(int_driver == 1) {
  75.         max_iter = int_nstep;
  76.         time = int_start;
  77.         time_step = (int_end - int_start) / int_nstep;
  78.     }
  79.  
  80.  
  81.     /* initialization of step for implicit symp algorithm */
  82.     /*
  83.     if(choice == 2) {
  84.         (int) f_p(vx1,1,vx,param,time,var_dim);
  85.         for(i=0; i<dim2; i++) delx[i] = time_step * vx1[2*i] + TINY;
  86.     }
  87.     */
  88.     
  89.  
  90.     /* main loop */
  91.     section_count = 0;
  92.     for (iter=1;iter<=max_iter && section_count < n_poincare;iter++){
  93.         /* to ensure the exactness of the ending time */
  94.         if(iter == max_iter -1){
  95.             if(int_driver == 1)
  96.                 time = int_end;
  97.         }
  98.         /* break if an orbit diverges */
  99.         for(i=0;i<var_dim;i++){
  100.             if(vx[i] > cutoff || vx[i] < -cutoff) {
  101.                 system_mess_proc(1,"Orbits appear to diverge off to an infinity! Stop!");
  102.                 goto done;
  103.             }
  104.         }
  105.         /* Not necessary since we are not using a newton's method
  106.         for iteration
  107.         if(choice==4){
  108.             for(i=0; i<dim2; i++){
  109.                 i2=i*2;
  110.                 vx1[i2] = vx[i2] + 0.9*delx[i];
  111.                 vx1[i2+1]=vx[i2+1];
  112.                 vx2[i2] = vx1[i2] + 0.2*delx[i];
  113.                 vx2[i2+1]=vx[i2+1];
  114.             }    
  115.         }
  116.         */
  117.  
  118.         int_onestep(vx1,vx2,vx,&time,time_step,var_dim,choice);
  119.  
  120.         /* from standard to window variables for plotting */
  121.         
  122.         /* Case: draw all orbits */
  123.         if(draw_all){
  124.             /* from window variables to full variables (full dimension)
  125.                include function variables if necessary */
  126.             to_full_variables(vfull,vx1,polar_coord);
  127.             /* test, plot and record an orbit */
  128.             (void) draw_record_orbit(vfull,iter,1);
  129.         }
  130.         /* Case: Draw Poincare sections */
  131.         else {
  132.             /* from window (or standard) to section variables
  133.                including function variables    */
  134.             to_full_variables(var_sec,vx1,polar_section);
  135.  
  136.             /* Test for intersections with the Poincare surface plane 
  137.                Surface is defined in terms of section variables
  138.                +function variables) */
  139.             if((var_sec_old[section_index] -section_constant) * (var_sec[section_index] - section_constant) <= 0 &&
  140.                 (int) (var_sec[section_index]-section_constant) == 0 ){
  141.                 if(linear_interpolation==1) {
  142.                     /* compute the time step from an old point
  143.                     to be on the section surface by linear
  144.                     interpolation */
  145.                     new_time_step = time_step * (section_constant - var_sec_old[section_index])/(var_sec[section_index]-var_sec_old[section_index]);
  146.                     for(i=0; i<var_dim; i++) vx2[i] = vx[i] + new_time_step / time_step * ( vx1[i] - vx[i]);
  147.                 }
  148.                 else if(linear_interpolation==2){
  149.                     new_time = time;
  150.                     new_time_step = time_step * (section_constant - var_sec_old[section_index])/(var_sec[section_index]-var_sec_old[section_index]);
  151.                     int_onestep(vx2,vx3,vx,&new_time,new_time_step,var_dim,choice);
  152.                 }
  153.                 else {
  154.                     /*
  155.                     system_mess_proc(1,"Error: Only linear interpolation available!");
  156.                     */
  157.  
  158.                     /* Only when Newton's method is used */
  159.                     new_time_step = time_step * (section_constant - var_sec_old[section_index])/(var_sec[section_index]-var_sec_old[section_index]);
  160.                     nl_interpol(&new_time_step,vx2,var_sec_old,var_sec,vx,time,time_step,var_dim);
  161.                 }
  162.                 section_count++;
  163.     
  164.                 /* from window variables to full window variables including function variables */
  165.                 to_full_variables(vfull,vx2,polar_coord);
  166.                 /* Test, draw, and record section */
  167.                 (void) draw_record_orbit(vfull,iter,1);
  168.             }
  169.         }
  170.         /* Reset new variables to old */
  171.         for(i=0; i<full_dim; i++) var_sec_old[i] = var_sec[i];
  172.         /* implicit symp algorithm ONLY */
  173.         if(choice == 4){
  174.             for(i=0; i<dim2; i++){
  175.                 i2=i*2;
  176.                 delx[i] = vx1[i2] - vx[i2];
  177.                 vx[i2] = vx1[i2];
  178.                 vx[i2+1] = vx1[i2+1];
  179.             }
  180.         }
  181.         else {
  182.             for(i=0; i<var_dim; i++) vx[i] = vx1[i];
  183.         }
  184.         /* break if error is detected or an interrupt is received */
  185.         if(stop) {
  186.             goto done;
  187.         }
  188.     }
  189.  
  190. done:
  191.     /* reset the final value of window variables */
  192.     to_window_variables(win_var_f,vx,polar_coord);
  193.  
  194.     if((int_driver==0 || int_driver==1) && int_algorithm==4){
  195.         sprintf(string,"Centered Euler: nok=%d nbad=%d\n",nok,nbad);
  196.         system_mess_proc(0,string);
  197.     }
  198.     else if(int_algorithm!=4 && linear_interpolation==0){
  199.         sprintf(string,"New Interpol: nok=%d nbad=%d\n",nok,nbad);
  200.         system_mess_proc(0,string);
  201.     }
  202.  
  203.     free_dvector(vx,0,var_dim-1);
  204.     free_dvector(vx1,0,var_dim-1);
  205.     free_dvector(vx2,0,var_dim-1);
  206.     free_dvector(vx3,0,var_dim-1);
  207.     free_dvector(var_sec,0,full_dim-1);
  208.     free_dvector(var_sec_old,0,full_dim-1);
  209.     free_dvector(vfull,0,full_dim-1);
  210.     free_dvector(delx,0,var_dim/2);
  211.     return(iter);
  212. }
  213.